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Abstract. Extending the ideal MHD stability code MISHKA, a new code, 
MISHKA-A, is developed to study the impact of pressure anisotropy on plasma 
stability. Based on full anisotropic equilibrium and geometry, the code can provide 
normal mode analysis with three fluid closure models: the single adiabatic model 
(SA), the double adiabatic model (CGL) and the incompressible model. A study 
on the plasma continuous spectrum shows that in low beta, large aspect ratio 
plasma, the main impact of anisotropy lies in the modification of the BAE gap 
and the sound frequency, if the q profile is conserved. The SA model preserves 
the BAE gap structure as ideal MHD, while in CGL the lowest frequency branch 
does not touch zero frequency at the resonant flux surface where m + nq = 0, 
inducing a gap at very low frequency. Also, the BAE gap frequency with bi- 
Maxwellian distribution in both model becomes higher if p± > py with a q 
profile dependency. As a benchmark of the code, we study the m/n = 1/1 
internal kink mode. Numerical calculation of the marginal stability boundary 
with bi-Maxwellian distribution shows a good agreement with the generalized 
incompressible Bussac criterion [A. B. Mikhailovskii, Sov. J. Plasma Phys 9, 190 
(1983)]: the mode is stabilized(destabilized) if py < p± (py > pj_). 


1. Introduction 

The magnetohydrodynamics (MHD) theory is widely applied in fusion plasma, 
providing a great aid in explaining various plasma instabilities and the plasma 
oscillating spectra below the ion cyclotron frequency. In modern toroidal magnetic 
confinement devices, the plasma contains significant fast populations originated from 
neutral beam injection (NBI) and ion cyclotron resonance heating (ICRH), inducing 
strong pressure anisotropy [I]. The magnitude of anisotropy can reach p\\ ^ l-7p± 
in a MAST beam heated discharge HE], or p± ss 2.5p|| in a JET ICRH discharge 
(4j, with p|| and p± the pressure parallel and perpendicular to the magnetic field 
lines, respectively. However, the physics of pressure anisotropy is not covered by the 
isotropic MHD theory. 

In the regime where wave-particle interaction is not important, a fluid approach is 
often used with a reasonable fluid closure (like the adiabatic condition for ideal MHD) 
for phenomena only related to the macroscopic quantities such as density, current 
and pressure. Many attempts have be made to incorporate anisotropy into the fluid 
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theory. Chew, Goldberger and Low (CGL) [5] first introduced the widely-applied form 
of pressure tensor and derived the double-adiabatic (CGL) closure, with its energy 
principle later derived by Bernstein et al [Bj. Unlike MHD, CGL assumes parallel and 
perpendicular pressures doing work independently in a collisionless plasma, therefore 
cannot reduce to MHD in the isotropic limit. It was found that CGL overestimates 
SW, the perturbed potential energy, compared to the kinetic theory, while MHD 
underestimates it 0EJ. Also, the mirror stability limit given by CGL does not 
match the result of kinetics theory mm- The major problem with CGL comes from 
the ignored heat flow when the mode frequency is comparable or smaller than the 
particle streaming frequency, especially in the vicinity of marginal stability boundary 
mm- Still, the CGL closure is implemented in many stability treatments, such as 
the ballooning modes ehee To overcome these drawbacks of CGL, some authors 
have proposed alternative fluid closures, for instance the double polytropic laws EE 
a higher-order-momentum closure EH EE and recently, the single adiabatic (SA) 
model [IS] which has the unique property of producing the same results as the MHD 
model for isotropic equilibria. Another pathway to overcome the drawbacks of CGL 
is to use hybrid approaches, in which thermal components are described by MHD 
and the fast ions by kinetics. The impact of pressure anisotropy is often investigated 
using kinetics energy princples mmm- In tokamaks, efforts have been made to 
study sawtooth modes (see Graves et al [22] and Chapman et al m and references 
therein) and interchange modes [22]. There are also significant developments in 
stellarators. The ANIMEC code [23] solves the 3D anisotropic equilibrium with the 
fast ion described by a guiding center distribution function, and is further applied to 
model anisotropy on LHP [24]. An energy principle which assumes non-interacting hot 
particles 23 is implemented in the ideal MHD code TERPSICHORE [2B] to model 
anisotropic-pressure interchange modes in a beam heated LHD discharge El!. Despite 
its shortcoming, the fluid approach can aid in the understanding of various effects due 
to its simple and intuitive nature. To date, there are few numerical studies on the 
oscillating spectrum of a toroidal anisotropic plasma. 

In the regime where significant wave-particle resonance exists, a pertubtive 
approach, in which the equilibrium and the linear mode eigenfunctions are modeled by 
fluid theory and the wave-particle interaction by kinetic theory, is widely implemented. 
In tokamaks, one of the most utilized tool chains is the HELENA-MISHKA-HAGIS 
combination E3 EH Eq], with the equilibrium, geometry and mode eigenfunctions 
calculated by ideal MHD, while the fast ion response and non-linear mode evolution 
are described by drift-kinetics equations. It has been successful in resolving the 
fast-particle-excited global Alfven eigenmodes (see reviews EU [32] and references 
therein). Recently, several equilibrium codes [3] [331 El] have been developed to 
study the equilibrium of anisotropic and toroidally rotating plasmas. For linear 
stability problem, efforts have been made to include the physics of diamagnetic 
drift and toroidal flow into MISHKA [33J EH for an isotropic equilibrium, while the 
impact of pressure anisotropy based on a full anisotropic equilibrium and geometry 
remains untouched. Our previous study using current remapping techniques shows 
that anisotropy can modify the q profile in MAST, inducing double TAE modes 
with different localization El EE and thus a double wave-particle resonance. This 
also serves as a motivation to develop a MISHKA-likc code to study the impact of 
anisotropy on linear stability, meanwhile drive a kinetic code using a fully anisotropic 
framework. 

This work is organized as follows. In Section^ we state our basic assumptions and 
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list the plasma equations used in the paper. Section [3] briefly describes the anisotropic 
equilibria and introduce the straight field line coordinates, serving as a base point 
for the stability treatment. Then in Section [|J we derive the linearized momentum 
equation, ideal Ohmic’s law and the fluid closure equations which are ready to use in a 
MISHKA-like numerical code. Section [5] introduces the implementation of the derived 
equations into a global normal mode code, MISHKA-A, and a continuous spectrum 
code, CSMISH-A. Using these tools, we study the impact of anisotropy on the plasma 
continuous spectrum and the internal kink mode, shown in Section [ 6 ] and Section 0 
respectively. We also compare the numerical results with existing analytical theory, 
serving as a code benchmark. Finally, Section [ 8 ] summarizes the paper and draws the 
conclusion. 

2. Plasma Model 

We start from a plasma described by the first two moments of the Vlasov Equation 
(the continuity and the momentum equation), the Maxwell Equations and the ideal 
Ohmic law. The basic equations are 


! + »(V.V) = 0, 

(i) 

p^ = -V-P+jxB, 

(2) 

f) R 

- = Vx(VxB), 

(3) 

j=Vxfl, 

(4) 

<1 

tn 

11 

0 

(5) 


where p is the mass density, V the mass velocity, P the second rank pressure tensor, 
j the current density and B the magnetic field. For simplicity, we use a natural 
MHD unit system where /To, the vacuum permeability, is set to 1. All electromagnetic 
fields, fluxes and vector potentials can be restored to S.I. units with a transformation 
■ /y/flo (e.g. B —X B/yfpf) and all currents with j —x y flloj- Equation dH) is 
the continuity equation. Equation ([2]) is the momentum equation. Equation ([3|, (J4]) 
and j5]) are the Maxwell Equations with ideal Ohmic law ignoring the displacement 
currents. The pressure tensor P takes the CGL form, i.e. 

p = p± 7 + abb, a = ^£ ( 6 ) 

with I the identity tensor, p± and p || the pressure perpendicular and parallel to the 
magnetic field, respectively. In our treatment, the finite Larmor radius (FLR) and the 
finite orbit width (FOW) effects are ignored. These effects can be important for fast 
particles, but resolving them requires FLR correction of non-diagonal pressure tensor 
terms (such as Chhajlani et al [3Sj for CGL) or kinetics/gyro-kinetics approaches, 
which are not considered in this paper. 

In this paper, we implement the standard linearization method, which expands all 
quantities into a combination of a time-averaging equilibrium part and a small time- 
dependent part, which varies with e . The mode frequency uj and growth rate 7 are 
related to A through the relationship A = 7 — ioo. By substituting these representatives 
into the plasma equations and considering the zeroth and the first order separately, 
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the equations are then converted into a time-independent equilibrium problem and a 
linearized stability problem. We drop the subscripts “0” for equilibrium quantities for 
convenience. 

To close the set of equations, one needs to introduce a “fluid closure” which relates 
PH and p± to other known variables. In this work, we examine three fluid closures: 
the single adiabatic model [18| , the double adiabatic model [5j, and the incompressible 
limit given by Mikhailovskii [39J . The single adiabatic model serves as a generalization 
of MHD. While keeping the adiabaticity assumption of MHD, it assumes that the 
parallel and perpendicular pressure are doing joint work, and therefore resolves the 
isotropic part of the pressure perturbation. This fluid closure equation is given by 


foil 

dt 


dp~± 

dt 


-v ■ v (|jP|| + fw.) 

- Qpu + ^p±) v • v - (Jp,i - b.(b. vn 


(7) 


in which the unit vector b = B/B is the direction of the magnetic field line. 
In contrast, the double adiabatic model assumes that parallel and perpendicular 
pressure do adiabatic work independently. The fluid closure equations, d/dt(p±/pB) = 
d/dt(p\\B 2 /p 3 ) = 0, after substituting Eq. Jl]) for dp/dt and B direction of Eq. ([3]) 
for dB/dt , are rewritten as 

^ = -v • Vp|| - P ||(V ■ V) — 2 p ,|b ■ (b ■ VV), (8) 

% = -V ■ V P± - 2p ± (V ■ V) + p±b ■ ( b ■ VV). (9) 

at 

Finally, the incompressible closure is obtained when the Lagrangian perturbed 
distribution function is set to zero, i.e. df /dt = df/dt+ V± ■ V/o = 0, where / 
is the Euler perturbed distribution function and fo is the equilibrium distribution 
function. After integrating over the velocity space, the incompressible closure fluid 
closure is given by 


dp~ll = _yl ( d P\\ \ 
dt \d S J B ’ 

ftpl = _T^1 ( dp± \ 
dt \ ds J B ’ 


( 10 ) 

( 11 ) 


where V 1 is the contravariant component of the straight field line coordinates (s, i9, <p), 
which will be introduced in the next section. 


3. Equilibrium and geometry 

For the zeroth order equilibrium problem, the time derivatives d/dt = 0. In this 
work, we ignore all equilibrium flows, i.e. Vo = 0. Using Eq. dSJ) in an axisymmetric 
tokamak geometry, the equilibrium magnetic field in cylindrical coordinate ( R , Z, <p) 
is written as 

B = V'k x Vp + FVtp, (12) 

where is the poloidal magnetic flux, F = RB V , and B v is the toroidal magnetic 
field. We note that unlike plasma with isotropic pressure, we do not require F to be 
a flux function. 
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Substituting Eq. ©, © and (1121) into Eq. m, the component in X7ip direction 
gives rise to a flux function F M f 3>) = RB V { 1 — A), while the V4/ direction gives the 
anisotropy modified Grad-Shafranov Equation (GSE). The modified GSE have two 
equivalent forms, the pressure form and the enthalpy form (See ;S] and references 
therein). In the pressure form of the GSE, the input profiles are specified by Fm(' h) 
and a 2D profile p\\(f&,B). This 2D pressure profile is usually obtained by taking 
the moments of guiding center distribution functions (401 either analytically (41) 
or numerically [42] (see Ref. [43] for a brief review). The enthalpy form of the 
GSE, when solved assuming the distribution functions are bi-Maxwellian and the 
parallel temperature is a flux function Ty = Tj|(\I/), requires four flux functions 
{H 1 Tj|, Fm, 0} as input, corresponding to the density, parallel temperature, toroidal 
field and anisotropy, respectively. The profile U(4') gives the radial shape of the 
density profile, and in isotropic plasma we have p = exp(iJ/Tj|). The profile 0(4') 
defines the anisotropy magnitude p±/p\\, which is given by 

P± B 


P\\ \B ~ 0Tj| | ’ 

The density and pressures are then linked to these profiles through 


P = 


P± H 

— exp—, 

Pll T ll 


and 


= pT\\, P± = pTj_ = pT r 


B 


B - 0Tii 


(13) 


(14) 


(15) 


These equation are identical to taking the moments of a bi-Maxwellian distribution 
function of the form in McClements et al 04] , written as 


F(p,E,*) = n r (*)- 


A{r) 


■ exp 


\E - pB 0 \ pB 0 


T,|(4/) 


(16) 


y / 2 *T ± (V) 

where A(r) is a normalization factor and 0 is just a convenient representation of the 
combination 

1 1 


0 UnW °‘ (17) 

In this paper, we will use this bi-Maxwellian model to explore the impact of anisotropy 
on stability, since it is the simplest model that captures pressure anisotropy for both 
ICRH and NBI. The model has limitations, such that it takes all species as a single 
bi-Maxwellian therefore cannot reproduce the long tail of ICRH fast ions, and that it 
omits any physics due to fine structure of pitch angle dependency of the distribution 
function (i.e. non-bi-Maxwellian structure). However, it does give the correct (py) 
and (p_ l), as well as A P\\/p\\, A p/p (the change of these profiles on a flux surface), and 
anisotropy A, which are not determined by a choice of the shape of the distribution 
function 00]. Here, (...) means flux surface average. We also mention that our stability 
treatment later on does not rely on the choice of equilibrium distribution function, as 
long as the modified GSE is solved self-consistently, and can provide T as a function 
of ( R , Z), i.e. the flux surfaces, for the stability treatment. 

The solution T(I?, Z) for the modified GSE is then mapped into the straight field 
line coordinates (s, $, <p), with s = \/4t and i? defined by 


d = 


B v dl 

qRBp 


qW = 


l 

27T 


B v dl 
RB n ' 


(18) 
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in which B p is the poloidal field and q the safety factor. The integrals are performed 
on a constant T surface clockwise facing the direction of e v and starting from Z = Z 0 , 
in which Z 0 is the Z coordinate for the magnetic axis. The metric coefficients of this 
curvilinear coordinate, gij and g lJ , as well as the Jacobian J, are defined by 


gij = Vx* • Vi 3 , 
J = yjdet(gij) = 


fqR 

F 


dr 

dx i 


dr 

dxi ’ 


(19) 

( 20 ) 


where / = d^/ds and det is the determinant operator, with(x, x , x 3 ) = (s,d,<p). In 
the straight field line coordinates, the contravariant equilibrium current is given by 


.! .2 

3 ~ J dti’ 3 J ds' 

.3 _ 1 ( d g 22 F d g 12 F\ 

3 J qR 2 dd qR 2 ) ’ 


( 21 ) 


B 1 = 0, B 2 = 


and the contravariant magnetic field components are given by 

F F 

qR 2 ' “ R 2 ' 

For the GSE with anisotropy in the straight field line coordinates, one can refer to 
Ref. ITS], as we will not restate it here. 


( 22 ) 


4. The perturbed equations in the straight held line coordinates 


In this section, we write our first order perturbed equations in the straight field 
line coordinates using contravariant and/or covariant representatives. Same as the 
original MISHKA, a set of “optimized” projections of V and B is used instead of 
the contra/co-variant projections. We use circumflexes to label these projections in 
order to distinguish them from the contra/co-variant projections, which are labeled 
by tildas. The perturbed fluid velocity V is expressed in its contravariant normal 
component V 1 , its binormal projection V 2 and its parallel projection K 3 , with 

V 2 = [VxB} 1 , = (23) 

The perturbed magnetic field B is calculated by taking the curl of the perturbed 
magnetic vector potential A (i.e. B = V x A). Then similarly, A is expressed in its 
covariant normal component Ai, its binormal projection A 2 and its parallel projection 
A 3 , with 


A 2 


[AxB ] 1 - A B 
B 2 ’ 3 - 


(24) 


The conversion between these projections and contra/covariant components of both 
V and A can be found in Ref. [22] and [15] . while B 1 , B 2 and B 3 are related to 
Ai, A 2 and A 3 through Eq. (90) to (92) in Ref. [18]. The covariant components are 
related to contravariant components through Bi = JfjgijBi. Finally, the perturbed 
magnetic field strength is given by B = B ■ b. 









Modeling the effect of anisotropic pressure on tokamak plasmas normal modes and continuum using fluid approaches7 
4.1. The ideal Ohm’s law 

Equation ([3]) , the ideal Ohm’s law, stays unchanged moving from isotropic plasma 
to anisotropic plasma. The equations are therefore identical to Ref. [23]: 


AAi = t> 2 , 

(25) 

AA2 = —Vi, 

(26) 

aa 3 = 0. 

(27) 


We recall that A = 7 — iui. When plasma equilibrium flow and resistivity are ignored, 
A 3 is an ignorable component, henceforth neglected. 


4-2. The momentum equation 


Perturbing Eq. ED, one obtains 


in which 


and 


dV — 

P0 - = ^- p + H , 


P = p~±{I-bb) +p\\bb + (p|| -p x ) I -fT b + b -ff 


H = (V x B) x B - B x (V x B). 


(28) 

(29) 

(30) 


The first two covariant components of H, Hi and Hi , are provided in Ref. [T5] and 
restated in | Appendix A| while H 3 is given in | Appendix A| as well. 

After some algebra, we reach the perturbed momentum equation covariantly in 
the straight field line coordinates: 


- B 1 B J 

XpVx = (1 - A)H 1 - d s p~± - dj{p\\ - p~± - 2ABB)-^- 
-A d s (BB) - (B j B 1 + B j BfjdjA 
-(pj| ~P~± ~ 2 ABB) -b + K^j , 

- Ik, IP 

A pV 2 = (1 - A )H 2 - dtip~± - dj(jp\ | - pi - 2A BB) ~ 
-Ad#(BB) - (B j B 2 + B j B 2 )djA 
~{p\\ ~ P~± ~ 2 ABB) • b + Ki'j , 


(31) 


(32) 


summing over index j = 1, 2,3, in which k = b ■ Vfo is the magnetic field line curvature 
with its covariant components and n 2 given in | Appendix A| Taking the dot product 
of Eq. (12811 with B, the third component of the momentum equation is written as 

Xp\B\ 2 V 3 = (1 - A )B j Hi - B J dj(p n - 2ABB) - AB j dj(BB) 

-djA(B j \B\ 2 + B j BB) - (pj, - pi - 2A BB){BV ■ 6), (33) 


summing over index j = 1,2,3. 
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f.3. The fluid closure equation 


For the single-adiabatic and double-adiabatic model, the fluid closure equations 
have similar forms in the straight field line coordinates, which are given by 

Apj, =-^[d 8 {JV 1 ) + d*{JV 2 )+d v (JV 3 )\-'y\\ 2 E 

-{v 1 d s + v 2 d»)f h (34) 

Api = -^[dsiJV 1 ) + d#{JV 2 ) + d^JV 3 )} - 7 ± 2 E 

-{v 1 d s + v 2 d»)f u (35) 

where 

n:i „ y2 

E=-d j (BV 3 )-V 1 Kl -—K 2 . (36) 

B fq 

For single-adiabatic model, we have 

14 2 2 

7||1 = 7-Ll = gP|| + jjP-L, 7||2 = 7-L2 = -P II - gP±, 

f\\ = f± = ^P|| + ^P±- (37) 

For double-adiabatic model, we have 

7||i=P|b 7||2 = 2p||, /||=P||, 

7-Ll = 2pj_, 7||2 = -p. l, f\\ = P. l- (38) 

There is no need to restate the incompressible fluid closure here, since Eq. m 
and m are already given in the straight field line coordinates. 

5. Numerical method 

Similar to the original MISHKA and its extension MISHKA-D/F, we use 
the following variables in our anisotropic extension of the MISHKA code, namely 
MISHKA-A (anisotropy): 

X\ = fqV 1 , X 2 =iV 2 , X 3 = iA u X 4 = fqA 2 , 

X 5 = ifv 3 , x 6 = fp~ ± , X 7 = fp\\. (39) 

These variables are then expanded poloidally and toroidally in Fourier harmonies with 
mode number m and n respectively, and radially in cubic/quadratic Hermite elements, 

i.e. 

oo N 

X a = e xt+inv ^2 (40) 

m——oo i/—l 

in which H u (s) is the cubic/quadratic Hermite elements and N the number of radial 
elements. The weak form is constructed by multiplying Eq. ©, © , ©, ©, 
©, © and © respectively by V l */{1 - A), V 2 */fq{ 1 - A), fV 3 */( 1 - A), A\/J, 
pq 2 A\/J, /p|| and fp±, converting the system into a linear algebra problem solving 

A Ni = Mi, (41) 


in which 


Ni=J2 I B{i,j)X*X j Jdsdd , 

3 =1 ’ 


( 42 ) 
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and 

8 /• BY* 

M r = J2J + A(i',j)M-i-X 3 

r) Y* rl Y 

+A(i,j')X*^. + A(i' ,j')—M-^-M-]Jdsdd. (43) 

We separate the matrix elements A(i,j) into 

A(i,j) = A°(i,j) + A A (i,j), (44) 

in which A°(i,j) are the common terms for MISHKA (isotropic) and MISHKA-A 
(anisotropic) and A A (i,j) are terms existing only in anisotropic plasmas. These matrix 
elements are given in |Appendix B| 

To obtain the continuous spectrum, we reduced MISHIKA-A to a continuum code 
(CSMISH-A). The method provided in Poedts et al [JS] (CSCAS) is implemented here, 
carrying the calculation in the vicinity of the singularity 4' —> Vto- 

6. Anisotropy impact on plasma continuous spectrum 

In this section, we study the continuous spectrum of an anisotropic plasma 
described by the SA and the CGL model, as well as the modification of anisotropy 
to the continuous spectrum. We present a set of examples with circular cross-section, 
large aspect ratio (e = 0.3) and low fj. The equilibrium solutions are computed 
by HELENA+ATF [3; using the enthalpy form of the modified GSE with the bi- 
Maxwellian distribution and the equilibrium thermal closure Tj| = Xj|(4/). We start 
from an isotropic MHD reference case with 

T(4' J v) = C 0 (1-4' J v) 2 + C 1 , RB v {* n ) = F q , 

H(* n )= < ^-( 1-* n ) 3 + C 2 , (45) 

where 'Pjv is the normalized flux surface defined as ’f'jv = 0 on axis and 'Pat = 1 at the 
edge, and Cq,Ci,C 2 and Fq are adjustable constants. Constant Fq indicates vacuum 
field strength. Constants C\ and C 2 are small values to make density and current 
profiles vanish at the plasma edge. The density and pressure profiles are given by Eq. 

(El) and (ED. The q profile monotonically increases from qo = 1.7 to <795 = 7. We 
choose ft = 1% on the magnetic axis. In the next step, we add anisotropy to this 
reference equilibrium. The 0 profile, which indicates the magnitude of anisotropy, is 
chosen to be constant. Therefore, anisotropy decreases from core to edge following 
the same trend of T, which is associated with on-axis beam heating or ICRH. For an 
individual anisotropic equilibrium, we specify a @ 0 , then iterate the Tj|, Fm and U 
profiles to keep (p*), (j) and (p) on each flux surface identical to the isotopic reference 
case. Here p* = (p|| + pj_)/2 and (...) means flux surface average. In this way the q 
profile and the metrics of these anisotropic equilibria are the same as the reference 
isotropic case to 0(e 2 (p|| — p_l)/p||). We have accordingly obtained equilibria ranging 
from p± = 1.7p|| (perpendicular beam or ICRH) to p || = 1.8p_L (parallel beam) at 
core. When we go to higher anisotropy like p± > 1.7p|| and py > 1.8pj_, we are unable 
to reduce the difference of qo between an anisotropic case (for example p± = 2 p\\) 
and its opposite case (p || = 2p±) to less than 1 % when we fix other parameters, 
since the flux surfaces of an anisotropic equilibrium is not completely reproducible 
by an isotropic equilibrium, or an anisotropic equilibrium with opposite magnitude of 







Modeling the effect of anisotropic pressure on tokamak plasmas normal modes and continuum using fluid approacheslO 


anisotropy [5]. Our start point is to identify the difference of anisotropic stability with 
equilibria in almost same conditions. Consequently, these higher anisotropy regimes 
are not explored here, because we are unable to keep them in these same conditions. 
However, our model and code are capable to describe cases with higher anisotropy, 
such as the pj_ = 2.5p|| discharge in JET. 

The continuous spectrum of these examples are then computed by CSMISH-A. 
Figure [Q shows the n = — 1 and m = 1,2,3 continuous spectrum of three cases : 
p± = 1.7p||, p_ l = p|| and p\\ = 1.8pj_ on axis, for (a) the SA model and (b) the 
CGL model. The linear growth rate of the continuous spectra in all these examples 
is observed to have 7 < lO -8 ^. We note that the small growth rate here is due to 
numerical errors (e.g. finite grid resolution) and is reduced by improving numerical 
precision. Therefore, we conclude that these continuous modes are stable. As in the 
ideal MHD spectrum, two sets of branches, a shear Alfvcn set (oj/ujao > 0.1) and 
a slow sound set (oj/loao < 0 . 1 ), appear at higher frequency and lower frequency, 
respectively. A resonance between m = 2 and m = 3 shear Alfvcn branches occurs at 
q = 2.5 surface and forms the TAE gap (Am = 1 gap) around s = 0.6. Meanwhile, 
a resonance between m = 1 and m = 3 forms the EAE gap (Am = 2 gap) at 
q = 2 surface around s = 0.4. The coupling between the shear Alfven and the slow 
branches forms the low frequency gaps (Am = 0 BAE gap). Moving to the edge, 
frequencies of the shear Alfven branches approach infinity as density approaches zero, 
while frequencies of the slow waves vanish as pressure goes to zero. 




6 

5 

4 

3 

2 

1 
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Figure 1. The n = —1, m = 1,2,3 continuous spectrum (left axis) and the 
q profile (right axis) of a plasma with (a) SA closure (b) CGL closure. The 
frequency is normalized to c^ao- the Alfven frequency on axis, and s = y/Wff is 
the standard flux label. The red dash line, black solid line and blue dash dot line 
shows respectively the cases with p± = 1.7p||, p 1 = p and p = 1.8px on axis. 
The EAE gap, TAE gap and BAE gap are labeled in each figure. 


Figure [T] also demonstrates the modification of anisotropy to the continuous 
spectrum. Anisotropy does not modify the main structure of the spectrum and the 
position of the gaps, but shifts the gaps and branches. For both models, around the 
core where the magnitude of anisotropy is higher, the difference between the three cases 
with different anisotropy is more significant. At the edge where anisotropy is vanishing, 
the three spectra merge to one. For the p\\ = 1.8pj_ case described by the SA model, 
all the shear Alfven branches are lowered (O.Olw^ on axis), while the slow branches 
are shifted up (7% on axis). For the CGL model, the lowest shear Alfven branch is 
almost unchanged, while the frequency of the slow branches increases by 14% on axis. 
The modification to slow branches will be investigated in Section 16.11 The change of 
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the shear Alfven branches can be explained by the change of these branches’ coupling 
to plasma compressibility through geodesic curvature, with different anisotropy and 
different model. Also, the q profile is only conserved to the reference isotropic case to 
0(e 2 (p || — p±)/p\\). With e = 0.3 in our example, the change of go is 0.01 (of 1.7) for 
the p\\ = 1.8p ± compared to the isotropic reference case, which will sightly modify all 
the branches. Looking at continuum gaps, the upper and lower accumulation points 
of both the TAE gap and the EAE gap are almost unchanged, meanwhile the upper 
accumulation point of the BAE gap is shifted up for both models (8% for SA and 4% 
for CGL). For the p± = l-7p\\ case, all the above modifications are reversed, with a 
similar magnitude of change. 

To understand the modification of anisotropy and the above differences, we study 
two specific feature of the continuous spectrum: its cylindrical limit and the low 
frequency BAE gap. The former one determines the main frequency of both the shear 
Alfven and the slow branches, and the latter describes the shear Alfven and slow 
coupling. 


6.1. the Cylindrical limit 


In the cylindrical limit, the equilibrium quantities are free of poloidal angle 
dependency. Therefore the coupling between two shear Alfven branches vanishes. 
Also, the geodesic curvature, which couples the shear Alfven branches and the slow 
branches, is zero. Building on Ref. |18l . we have computed the continuum in the 
cylindrical limit. We retain the ignored (pu — pif){b\b + bbf) term in the perturbed 
pressure tensor in Ref. [18] , therefore the missing firehose factor 1 — A for the single- 
adiabatic Alfven branches is now recovered. The frequency of mode (m, n ) is now 
simply given by 


2 _ 2 _ (1-A )B\ , 

w A,SA — w A,CGL — 7712 ( m /9 ' 


pR i 


■nf 


^S,SA — T 


+ b± 


0 

B 2 


%S,CGL ~ 


3P11 + 3P-L + R 2 pRo 
3p|| B 2 


(m/q + n) 2 


(m/q + n) 


(46) 


(47) 


(48) 


2 p± + B 2 pR 2 

where Rq is the major radius of the magnetic axis. Here, “A” in the subscript labels 
the shear Alfven branches and “S” labels the slow branches. Inspection of Eq. 051) 
shows that the cylindrical shear Alfven continuum is not fluid closure dependent. The 
anisotropy modifies these branches by the firehose factor 1 — A. This is consistent 
with previous results [106]. In contrast, the slow branches, as shown by FigQ] have 
strong fluid closure dependency and anisotropy dependency, with u>s,sa (jJS,CGL 
even when the equilibrium is isotropic. In the isotropic limit, the SA model reduces 
to the result given by ideal MHD with adiabatic gas law, while the CGL model does 
not converge to ideal MHD. Indeed, the frequency ojs.CGL is roughly 35% larger than 
us,sa when the plasma is isotropic. As in Eq. (071) and 051) . the frequency of the 
slow branches with both model are increasing when p\\/p± increases, if (p*) is kept 
constant, although CGL model shows more significant change compared to SA. We 
have compared the result from CSMISH-A in the cylindrical limit (very large aspect 
ratio) with Eq. 06l) to 051) for both SA and CGL, showing very good agreement. 
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6.2. The BAE gap change due to anisotropy 

The low frequency gap (BAE gap) 07] appears on the resonant flux surface where 
m + nq = 0, and is induced by the finite compressibility of the plasma. Inspection 
of FigJT] shows that for different magnitude of anisotropy, the width of this gap is 
changed. Also, the gap width is different for the SA and the CGL model, implying 
its dependency on fluid closure model. Figure [2] zooms in into the q = 2 BAE gap 
in FigJT] for the anisotropic case with p\\ = 1.8p_L on axis. Only the major m = 2 
shear Alfven branch and the m = 2 ± 1 slow side bands are shown here. In Fig[7] 
the frequency of the upper, middle and lowest branches on the resonant flux surface 
(located at s = 0.38) are labeled as W 3 , w 2 and uq respectively. The BAE gap of the 
SA model has the same structure as an isotropic plasma described by the MHD model. 

Its lowest branch approaches zero when m + nq = 0, i.e. ui 1 = 0. To the contrary, in 
CGL we have uq > 0, inducing an additional gap at very low frequency. 




Figure 2. Zooming into the q = 2 BAE gap of the n = — 1 continuous spectrum 
of the anisotropic case in Fig^with p n = 1.8pj_ on axis, for (a) the SA model 
and (b) the CGL model. The blue dash lines are the incompressible m = 2 shear 
Alfven branch. The vertical lines indicate the flux surface where q = 2, and the 
incompressible m = 2 shear Alfven branch hits zero. The red solid lines are the 
coupled m = 2 shear Alfven branch and m — 1,3 slow branches due to finite 
compressibility (SA or CGL), with the frequency of the upper, middle and lowest 
branches labeled as oj 3 , LJ 2 and uji at q = 2 surface, respectively. 


In this section, we are only interested in 0 ^ 3 , the upper accumulation point of a 
BAE gap, which determines the gap width. We study two separate cases, with the 
gap located at a low q position (q = 1.33) and a high q position (q = 3), as shown in 
FigJ3](a) and (b), respectively. The frequencies in Fig J3] are normalized to the analytic 
ideal MHD value of W 3 for the reference isotropic case [35], written as 

= (7 P + B*)pRi ( X + 2 ?) ’ (49) 

with 7 = 5/3. Figure [3] (a) shows that for q = 1.33, the SA closure gives a greater 
ui 3 when py > p±, and a smaller ui 3 when py < p±. It’s almost a linear function 
of (p || — p±)/p*. The change of W 3 is roughly 8 % for p || r; 1.5pj_ or p± « 1.5p||, 
the farthest right and left data points in the figure. For the CGL closure, W 3 is 7% 
higher than the isotropic ideal MHD reference case. It’s dependency on (py —pi_)/p* 
is almost negligible. Moving to Fig[3](b) where q = 3, in SA model the dependency of 
to 3 on (p || — pj_)/p* becomes higher, with a 12% change for p\\ ki 1.5pj_ or p± « l-5p||. 
Meanwhile, the ratio ojs/oj^mhd decreases to 1.03 in the isotropic case, and the UJ 3 
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for CGL has a weak dependency on anisotropy: about a further 5% change for the 
extreme cases. 
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Figure 3. The change of the BAE gap upper accumulation point frequency (^ 3 ) 
due to the change of anisotropy for a BAE gap with (a) q = 1.33, n = —3 (b) 
q = 3.00, n = — 1. The local magnitude of anisotropy is described by the relative 
difference of py and p±, i.e. (py — p±)/p*. The frequency of CJ 3 is normalized to 
the analytic ideal MHD value of u )3 for the reference isotropic case, as shown by 
the horizontal dash line. The symbols are numerical results from the CSMISH-A 
code: blue squares and solid lines for SA, red circles and solid lines for CGL. 


7. Anisotropy impact on the internal kink mode 


In this section, we study the impact of anisotropy on the n = 1 internal kink 
mode in a tokamak plasma with large aspect ratio (e = 0.1) and circular cross section. 
This also serves as a benchmark of MISHKA-A working as a global normal mode code. 
For simplicity, the equilibrium distribution function is taken to be bi-Maxwellian. 

We start from a reference isotropic equilibrium with the current profile and the 
pressure profile taking the form, 


O') = Jo(l - y N ), (50) 

P\\ =P± =Po(l - ^n), (51) 


where jo and po are constants. The density profile is taken to be constant, i.e. p = po- 
The safety factor on axis, qo , and the ratio of kinetic energy to magnetic energy, /3, can 
then be adjusted by changing the ratio po/jo and the vacuum field. The safety factor 
q is monotonically increasing: only one q = 1 surface exists in the plasma. Similar to 
Section O based on this reference isotropic case we change the 0 profile with 0 = 0 O 
in our equilibrium code HELENA+ATF, meanwhile keeping (p*) = ((p±) + (j?||))/2, 
0) and {p} unchanged. In such a way the q profile and metrics are identical to our 
reference isotropic case to 0(e 2 ). The relative anisotropic profile is then approximately 
given by 


(P- l) = 1 

(P||) 1 0.(1 Tjv) 1 


(52) 


with which the magnitude of anisotropy peaks on axis and vanishes at the boundary. 
Here a is an adjustable constant proportional to 0o- 

In the incompressible limit, the plasma kinetic response to the perturbation is 
ignored. The stability of the internal kink mode is determined by the sign of the 
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perturbed fluid toroidal potential energy SWt■ When SWr < 0, the plasma is 
unstable. According to the analytical calculation of Bussac et al 09] and Mikhailovskii 
[391151)] . the stability criterion of the n = 1 internal kink in such a scenario, namely 
the generalized Bussac criterion, is described by 


Sw + P pA > 0, 


(53) 


where Sw is a quadratic function of the value of /3 P on the q = 1 surface, with the 
coefficients determined by the q profile. The quantity Bussac /3 P , as a indication of 
the pressure gradient, is defined as 


= 


2[p(T) -p(T)] 

B 2 P m 


(54) 


where p is the average pressure inside the certain flux surface, i.e. 


p{^ l) = f pdS/ f dS. 


(55) 


For anisotropic plasma, /3 P is replaced by /3 P * = (/3 p y + /3 p± )/2. The second term in 
Eq. m . fd p A: is obtained from Eq. (1591) replacing p by (p|| + p± + c)/2, and taking 
the value on the q — 1 surface as well, where c is defined through partial derivative of 
p± as 


B 


(w) - 


= 2pj_ + c. 


For a bi-Maxwellian plasma, c is simplified to 

2 p± 2 


CbM = — - 


P II 


(56) 


(57) 


The generalized Bussac criterion takes into account only the lowest order of 
the poloidal variation of py and p~± , and neglects the shaping effect EP of pressure 
anisotropy, leading to its discrepancy from full numerical results when the fast particle 
distribution function has strong and/or complicated poloidal dependency (e.g. with 
neutral beam heating) [52]. The bi-Maxwellian plasma we use here only has a weak 
poloidal dependency, satisfying the use of the generalized Bussac criterion. From Eq. 
(1571) . /3 p a is positive when py > p± and negative when py < p±. We would expect 
the plasma to become less stable compared to the reference isotropic case if p± > py 
(a > 0 ) and more stable when p± < py (a < 0 ). 

To obtain the marginal stability boundary numerically, we plot the internal kink 
growth rate as a function of /3* for different a in Fig Q] (a). Figure[4] (a) shows that 
in anisotropic plasma, same as Bussac et al , the linear growth rate of the internal 
kink mode increases with /3*. For the same /?*, the growth rate is higher when a 
becomes more positive. On the other hand, the growth rate is reduced, or the mode is 
stabilized, when a becomes more negative. This is in agreement with the prediction 
of the generalized Bussac criterion. 

The critical /3* at marginal stability is extrapolated from Fig[4] (a) by fitting 7 
into a quadratic function of /?* and obtaining the fitted curve’s intersection with the 
x axis. Picking different qo and different a, the marginal stability boundary is then 
plotted in Fig[4](b) with a comparison against Eq. (1531) . Figure[4](b) shows that when 
a = 0, i.e. the plasma is isotropic, the stability limit given by MISHKA-A is in good 
agreement of the analytical Bussac limit. When a > 0 (pj_ > py), the anisotropic 
incompressible fluid force is destabilizing, reducing the required pressure gradient to 
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drive the instability. On the other hand, if a < 0 ( p± < py), the anisotropic geometry 
is stabilizing. We note that when q 0 is close to unity, the stabilizing/destabilizing effect 
is greater, pushing the stability limit further from the original Bussac limit. This is 
due to the fact that when qo is close to 1, the first term in Eq. (1531) . Sw , is smaller. 
Therefore a tiny change in /3 p a will lead to a dramatic impact of the stability limit. 
We also note that the magnitude of anisotropy in Figj4] is small (with pj_ = 1.25py on 
axis for a = 0.2, or py = 1.2pj_ on axis for a = —0.2). We would thus expect that a 
moderate or large anisotropy will have a much greater impact to the n = 1 internal 
kink mode. 

We observe a small discrepancy between the generalized Bussac criterion (lines) 
and the numerical result (symbols) in Fig|4|b) for the a < 0 cases. One possible 
reason is that in the derivation of the generalized Bussac criterion, the eigenfunction 
is assumed to stay the same as the isotropic reference case. Also, the perturbed 
parallel electric field B and the perturbed parallel flow V ■ b are ignored. These 
neglected features, when taken into account numerically, may have some impact on 
the marginal stability limit. Nevertheless, Fig|4] (b) gives a fairly good benchmark of 
the MISHKA-A code. 




Figure 4. (a) The growth rate of the n = 1 internal kink mode as a function 
of ft* 2 for a plasma with qo = 0.9. The parameter a determines the magnitude 
of anisotropy, with pj_ > py for a > 0 and pj_ < Py for a < 0. The growth rate 
7 is normalized to Alfven velocity Va • (b) The modified Bussac critical ft* as a 
function of qo for different anisotropy magnitute a. The lines are analytical result 
calculated from Eq. lf53l l and the symbols are numerical results extrapolated from 
(a). 


The above treatment ignores the compressional response of the plasma and keeps 
only the incompressible part. According to the kinetic theory, the compressional 
response can either be stabilizing or destabilizing, depending on the fast particle 
distribution function, the diamagnetic effects, FLR/FOW effects and other non-ideal 
effects (see for example the review of Graves et al [20] and Chapman et al [21]). A 
full treatment of the n = 1 internal kink mode will require a 5f method and possibly 
the involvement of a kinetic code. Nevertheless, we can still conclude on that the 
anisotropic incompressible fluid force of a plasma with p± < py (a < 0) is more stable 
than its isotropic counterpart, and therefore needs less stabilizing effects from kinetic 
response to stabilize, while a plasma with p± > py (a > 0) needs more. 

Finally, we investigate the compressional response of a plasma described by the 
CGL model. We couldn’t find any unstable modes for our choice of current and 
pressure profile, despite a scan across parameters 0.6 < qo < 1, 0 < (5* < 0.5 and 
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0.5 < pj_/p\\ < 2. It’s long been known that for isotropic plasma we have [UH] 

SWmhd < SWk < SWcgl , (58) 

where 5Wk is the perturbed potential energy given by the kinetic theory. For 
anisotropic plasma, although not rigorously proved, it is very likely to have SWk < 
SWcgl ■ With the CGL gives a prediction that the plasma is stable, we can conclude 
that for our choice of profiles and parameter space, it is possible to stabilize the internal 
kink mode by plasma compressional response. 

8. Conclusion 

We derived and implemented the linearized fluid equations with anisotropy in 
the straight field line coordinates based on three fluid closures: the double-adiabatic 
model (CGL), the single-adiabatic (SA) model, and the incompressible model. The 
ideal MHD normal mode code MISHKA has then been extended to its anisotropic 
pressure version, MISHKA-A (and the continuous spectrum code, CSMISH-A). Using 
these numerical tools, we find that anisotropy mainly modifies the continuous spectrum 
by changing the slow branches and the BAE gap. The change of the slow branches 
is in accordance with the analytical result, with a different prediction for the SA 
model and the CGL model. For the BAE gap, the lowest branch touches zero at the 
resonance flux surface for SA/MHD, but does not for CGL. Meanwhile the change in 
frequency of the upper accumulation point depends on the local q value, the magnitude 
of anisotropy and the fluid closure. Finally, we study the impact of anisotropy to the 
internal kink mode numerically. If only the incompressible fluid force is considered, 
we find that for a bi-Maxwellian plasma, the marginal stability boundary is in good 
agreement with the analytical result of Bussac et al and Mikhailovskii: the plasma 
is stabilized if p± < p\\ and destabilized if p\\ > p±. Also, a parameter scan reveals 
that for our choice of profiles the internal kink mode is stable, if the CGL closure 
is implemented. This indicates the possibility for these modes to get stabilized by 
the plasma compressional response, and that CGL is too strong for the estimation of 
instabilities. 

In this work we restrict our study to large aspect ratio, low beta plasma, 
when the equilibrium can be reproduced similarly by an isotopic equilibrium with 
an 0(e 2 (p|| — p±)/p\\) difference. In the future, we plan to study the impact of 
anisotropy on global eigenmodes, and the possibility of using these eigenmodes as 
MHD spectroscopy to infer pressure anisotropy. For example, as indicated by the 
change of the BAE gap due to anisotropy, the corresponding modification to a 
global BAE may serve as an estimation of pressure anisotropy or a validation of 
the fluid closure model. We also plan to investigate tokamak plasmas with high 
/3, low aspect ratio and large anisotropy, where the current profile and q profile 
are dramatically modified by anisotropy, and where the anisotropy shaping effect 
is important. Finally, we plan to study experimental data from MAST, with the 
equilibria reconstructed anisotropicly by the EFIT-TENSOR code [Mj, and compute 
the wave-particle interaction. 
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Appendix A. Auxiliary formulas 

Here we present the formular for Hi, H 2 , H 3 and covariant components of the 
magnetic field line curvature n: 

Hi = J{j 2 B 3 - j 3 B 2 ) - J^dsiguB 1 + g 22 B 2 ) 

+ + qd^ignB 1 + g^B 2 ) - -^d s (R 2 B 3 ), (A.l) 

H 2 = JifB 1 ~ j'B 3 ) + -^d v (g v2 B l + g 22 B 2 ) 

-JL ds {R 2 &), (A.2) 

H 3 = J(fB 2 - fB 1 ) - B 2 dfg 12 B 1 + g 22 B 2 ) 

+B 2 B 3 8 dR 2 , (A.3) 

F / d q\\/^\ 2 dF <9 V'E-VjA 
Kl ~ ~ qBR 2 Vds BF + Q dsB + BF ) ’ 

F d (F\ 

K2 ~ ^~RfB~dd \73 ) ’ 

«3 = -• 

q 

Appendix B. Matrix elements 

Appendix B.l. The momentum equation 

The left-hand sides matrix elements 13(1,1), 73(1,2), 73(2,1) and 73(2,2) are 
identical to those given in the appendix of Ref. [35] dividing by 1 — A. Elements 
73(1,5), 73(2, 5) and 73(5, 5) are given by 


73(1,5) = ipo-$rj=, —-V??, 
r £m 

(B.l) 

B(2,5) = Po ^-\V^\ 2 , 

(B.2) 

B(5,5) = Po ^-\B\ 2 . 

(B.3) 


(A.4) 

(A.5) 

(A.6) 
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For Eq. (1311) . the matrix elements A°(l, 3), A°( 1', 3), A°(l, 4), A°(l,4'), A°( V, 4) 
and A°(l',4') are same as those in the appendix of Ref. [33 and [3B], except that 
dF/ds in Ref. 33 and ES] is now replaced by dF/ds. Other A°(i,j) elements coming 
from Eq. (ED are 


R 2 


W-V-Tf? 

= ;!(£) 7-• 4 ” (1 - 7) ' 


d BR 2 d V'F- W 

* *dd F\B\ 2 


+i(m + nq) 


V'h ■ S7d 
F m \B\ 


2 ’ 


(B.4) 

(B.5) 

(B. 6 ) 


For Eq. (13^1) . the term A°(2,4) is same as Ref. [33 > but again changed its dF/ds 


terms to dF/ds. Other elements are given by 

. n , 1 , „ 0 o 9 , (rfi — m)m „ 

A°(2, 3) = - —— (mmF 2 + nV|V>F| 2 ) - - - —-— F, 

JQf jq 

(B.7) 

A°(2,4') = -i-(mF 2 - ng|W| 2 ) + 

fqF fq 

(B. 8 ) 

A°(2,6) = ^--A°(2,7), 

JrM 

(B.9) 

A°(2,7) = * (|Vtt| 2 d#B - F 2 d#B + FBd$F) 

J r m\B\° 



(B.10) 


. . |V ^| 2 

+,m + nq) 7^W 

Also, A°(i,j) elements from right-hand side of Eq. (1331) are listed as following : 
F OF 


A 0 (5,3) = i(m + nq) 
A 0 (5,4) = 


qR 2 dd ’ 
m + nq d\ V\I /| 2 m + nq 


(B.ll) 


qR 2 ds 
+ (m + nq) 


q 2 R 2 F 
fF 3 Vf- Vd 
qR 2 dd F 


dq OF 


l y ^| 2 (F-A-g 


ds 


+ (m + nq) 


F dF 
qR 2 ds 


—I 


A 0 (5, 6 ) = -i 


F dF dq 
q 2 R 2 dd ds ’ 
1 dB 


(1 - A )B dd 
0 / 


^°(5,7) = ^±^-A°(6,5). 

The anisotropy related terms are given by 
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(B.12) 

(B.13) 

(B.14) 
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A a (1A) = 


dp A 
Fm 
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and 
hi = - 

hi = 


F 2 


fqR 2 ’ 
VT ■ Vd 
qR 2 : 


ho = 


|VT| 2 


hr. = 


fR 2 
|5| 2 


h% = - 


|W | 2 rfg 

fq 2 R 2 ds' 


fq 


(B.30) 


Appendix B.2. The ideal Ohm’s law 

For the ideal Ohm’s law equations (Eq. C5l) and (I2S1) ), we have 
B(3,3) = ^°(3,2) = l, 

5(4,4) = -A°(4,l) = l, 

Appendix B.3. The single/double- adiabatic fluid closure equations 


(B.31) 

(B.32) 


The matrix element B(6 ,6) and 5(7, 7) are identical to Ref. j[5([! 5(7, 7). For the 
single/double-adiabatic model Eq. (l34l) and (l35l) . the A°(i,j) elements are given by 
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(B.34) 


+°(6,2) =7_L 1T ^-(n^^- m ) +*TT7? vrjC/x -7 xi) 


\B\ 2 

+*7X2 


F 2 
1 d F 
BddB' 


\B \ 2 dd 


d 


^4°(6, 5) = —7_Li(m + nq) - 7_l 2 /n + nq) + *-^(/j_ - 7 x 1 ) 


dd 


1 dB 
+n±2 B ~dd’ 


(B.35) 


(B.36) 


Replacing f± by /||, 7 _li by 7 p and qj _2 by 7 || 2 , we will reach the matrix elements 
A 0 (7, 1'), +°(7,1), A°{ 7,2) and +°(7,5). 


Appendix B.f. The incompressible fluid closure 

The matrix element 5(6,6) and 5(7,7) are identical to Ref. 
A°(i,j) elements originated from Eq. (flUl) and (fill) are given by 


+ u (4,l) = - — 


R 2 f dp± d$p± dB 




F \ ds d$B ds 
R 2 f dp\\ p\\ - p± dB 
B ds 


5(7,7). The 

(B.37) 

(B.38) 
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